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1.  Introduction  and  Background 


The  primary  purpose  of  anti-vehicular  landmines  is  to  inflict  damage  on  an  Army  vehicle  and  its 
occupants.  Since  the  blast  originates  from  the  ground,  the  occupants  experience  the  full  effects  of 
the  blast  through  their  lower  legs,  the  part  of  the  body  closest  to  the  mine.  By  studying  the  effects 
of  anti-vehicular  landmines  on  the  human  body,  modifications  may  be  made  to  the  vehicles  to 
mitigate  bodily  harm. 

The  closest  simulation  of  an  anti-vehicular  landmine  detonation  is  an  experiment  that  uses 
explosives  and  human  surrogates.  The  entire  setup,  including  the  dummy,  sensors,  and  test  stand, 
is  expensive.  The  experiment  could  take  weeks  to  prepare,  the  sensors  could  fail  to  read  the  data, 
and  the  dummy,  along  with  everything  else,  could  be  destroyed  in  a  single  test  without  collecting 
useful  data.  In  this  sense,  it  is  logical  to  create  a  computer  model  of  the  events  of  a  landmine 
event. 

The  first  step  in  simulating  a  landmine  event  is  to  use  a  finite  element  model  to  find  a  solution 
numerically.  In  this  instance,  the  solution  would  be  the  force  experienced  by  the  tibia  in  the 
vertical  direction.  However,  since  a  model  of  the  event  would  require  significant  computational 
power,  a  finite  element  model  would  not  be  ideal  for  a  quick  approximation  of  the  problem. 

The  ultimate  objective  is  to  create  a  simplified  computer  model  that  would  be  cost  efficient  in 
approximating  the  force  experienced  in  the  lower  leg,  which  could  then  be  used  to  detennine  the 
resulting  leg  damage.  This  report  describes  in  detail  the  process  in  which  the  simplified  computer 
model  was  developed  and  its  accuracy  in  predicting  the  forces  experienced  by  the  lower  leg  in  an 
anti-vehicular  landmine  blast  was  verified. 


2.  Experiment  and  Calculations 


The  goal  is  to  relate  the  floor  displacement  to  lower  leg  force  and  verify  the  analytical  solution 
using  the  results  from  testing  the  Hybrid-Ill  dummy  in  the  Test  Rig  for  Occupant  Safety  Systems 
(TROSS)  system  presented  in  Bir  et  al.  (2006).  Given  a  floor  displacement  over  a  period  of  time, 
an  analytical  model  should  be  able  to  accurately  predict  the  maximum  force  experienced  by  the 
tibia  of  the  occupant. 

2,1  Model 

This  analytical  model  uses  a  series  of  springs  and  masses  to  model  the  tibia.  Beginning  with  a 
one-spring,  one-mass  system,  the  model  has  been  expanded  to  calculate  the  maximum  tibia  force 
for  a  given  number  of  springs  and  masses  (figure  1). 
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Figure  1.  Analytical  model,  showing  a  (a)  one-mass,  one-spring 
system  and  (b)  an  n  mass,  n  spring  system. 

For  the  model,  the  n  masses  are  distributed  at  constant  intervals  and  labeled  from  1  to  n  starting 
at  the  knee  and  ending  at  the  ankle.  The  bottom  of  the  boot  is  rigidly  attached  to  the  floor.  All  n 
masses  have  one  degree  of  freedom  in  the  vertical  direction.  The  springs  are  assumed  to  be  ideal 
and  follow  Hooke’s  Law. 

This  analytical  model  uses  the  SI  unit  system. 

2.2  Calculations 

For  each  mass,  the  vertical  coordinate  is  calculated  using  a  combination  of  basic  physical 
concepts  and  a  numerical  integration  method.  The  masses  are  calculated  from  the  geometric 
parameters  provided  by  the  Hybrid-Ill  LS-DYNA  model.  Since  there  are  a  finite  number  of 
masses  in  the  model,  the  tibia  mass  is  discretized  into  n  parts  and  distributed  unifonnly. 

The  initial  conditions  for  the  vertical  coordinate,  vertical  velocity,  and  vertical  acceleration  are 
initialized  to  assumed  values.  For  this  model,  the  vertical  coordinates  are  initially  distributed 
evenly  over  the  length  of  the  tibia  given  the  number  of  masses  desired.  Two  masses  are  placed 
below  the  tibia,  one  at  the  ankle  and  the  other  at  the  heel.  Both  the  vertical  velocity  and  the 
vertical  acceleration  are  initialized  at  0  m/s  and  0  m/s  ,  respectively.  The  decision  to  leave 
gravity  out  of  the  problem  is  based  on  the  fact  that  the  accelerations  the  masses  experience  are  an 
order  of  magnitude  greater  than  that  of  gravity. 
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Using  the  forward  Euler  method,  the  midpoint  method,  and  the  two-step  Runge-Kutta  method  for 
numerical  integration,  the  values  for  the  states  for  the  previous  time  step  were  used  to  extrapolate 
the  value  for  the  states  in  the  current  time  step.  In  order  to  calculate  the  vertical  acceleration  for 
each  of  the  masses,  a  basic  force  diagram  was  used,  resulting  in  the  following  equations: 

Node  1 :  Knee:  mi  yf’  =  k(fi  -  (yi  -  y2))  ( 1) 

Nodes  2  to  n-1 :  Intermediate:  m,  y,”  =  k(l,  -  (y,  -  y,+,))  -  k(li_i  -  (y,_i  -  y,))  (2) 

Node  n:  Ankle:  mn  yn”  =  k(ln  -  (yn  -  yfl00r))  -  k(ln_i  -  (y„_i  -  yn))  (3) 

where  1;=  original  length  of  spring  i,  and  y;  =  vertical  coordinate  of  mass  i. 

Using  Hooke’s  Law  and  the  defonnation  of  the  tibia  in  the  vertical  direction,  the  total  force  on 
the  tibia  for  each  time  step  is  calculated. 

Since  the  model  does  not  implement  dampers,  the  results  from  this  analysis  oscillate.  The 
maximum  tibia  force  value  is  the  first  local  maximum  of  the  force  vs.  time  graph. 

2.3  MATLAB  Script 

Given  a  *.dat  file  for  the  floor  displacement  over  a  certain  period  of  time,  MATLAB  loads  the 
data  into  a  matrix.  Lrom  this  file,  the  time  and  corresponding  floor  displacement  are  stored  in 
column  arrays.  The  number  of  masses  is  stored  in  the  variable  n,  and  the  masses  and  lengths  of 
the  tibia  and  femur  are  stored  in  variables  as  well. 

Using  three  matrices,  the  state  variables,  vertical  coordinate,  vertical  velocity,  and  vertical 
acceleration  are  stored  for  each  mass  and  each  time  step.  An  iterative  loop  is  used  to  populate  the 
matrices  with  the  appropriate  value  calculated  via  the  method  described  in  section  2.2. 


3.  Results  and  Discussion 


In  order  to  verify  that  the  analytical  model  created  in  MATLAB,  two  sets  of  simulated  data  were 
produced  from  an  LS-DYNA  simulation  of  a  detonation  occurring  beneath  a  plate.  The  simulated 
results  are  compared  to  the  experimental  data  provided  in  Bir  et  al.  (2006).  Running  the  script 
using  the  original  data  with  two  springs  led  to  the  results  listed  in  table  1 . 

Table  1.  Values  calculated  from  simulated  data  using  analytical  two-spring  model. 


Data  File 

Max  Displacement 
(m) 

Max  Velocity 
(m/s) 

Max  Tibia  Force 
(N) 

Disp2.dat 

0.7910 

244.3158 

3,081,000 

Disp3.dat 

0.3555 

45.8794 

371,100 
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Since  both  of  the  simulated  data  sets  provided  maximum  floor  displacements  and  maximum 
velocities  that  do  not  lie  within  the  range  of  the  experimental  values,  the  data  sets  were  scaled  to 
provide  a  maximum  floor  displacement  within  the  range  of  the  experimental  data.  The  script  was 
run  again  with  two  springs  and  a  scaling  factor  for  each  of  the  two  data  sets  (table  2). 

Table  2.  Values  calculated  from  the  scaled  simulated  data  using  analytical  two-spring  model. 


Data  File 

Scale 

Max  Displacement 
(m) 

Max  Velocity 
(m/s) 

Max  Tibia  Force 
(N) 

Disp2.dat 

1/46 

0.0172 

5.3112 

66,970 

Disp3.dat 

1/21 

0.0169 

2.1847 

17,670 

From  this  analysis  using  the  forward  Euler  method  of  numerical  integration,  it  may  be  noted  that 
the  scaled  data  from  DISP3,  in  terms  of  maximum  displacement  and  maximum  velocity,  most 
accurately  matches  the  experimental  data.  For  the  rest  of  the  analysis,  only  1/21 -scaled  DISP3 
data  were  used  for  the  purposes  of  verification.  Despite  similar  values  for  displacement  and 
velocity,  the  maximum  tibia  force  was  over  two  times  that  of  the  experimental  data. 

In  the  second  round  of  analysis,  using  only  the  1/21-scaled  DISP3  data,  the  tibia  force  was 
computed  using  2  to  8  springs  with  distributed  masses  along  the  length  of  the  tibia.  Figure  2 
shows  the  tibia  force  versus  time  for  the  9-spring,  forward  Euler  analytical  model. 


Figure  2.  Tibia  force  versus  time  for  the  9-spring,  forward  Euler  analytical  model. 
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As  the  tibia  is  further  divided  into  more  springs  and  masses,  the  oscillatory  behavior  masks  the 
first  maximum  tibia  force,  so  the  maximum  number  of  springs  and  masses  that  may  be  used  in 
the  analytical  model  is  8. 

The  solution  to  the  n-spring  model  is  then  calculated  using  the  midpoint  method  for  numerical 
integration.  The  oscillatory  behavior  that  appeared  in  the  9-spring  forward  Euler  model  appeared 
in  the  27-spring  model  using  the  midpoint  method,  allowing  for  the  tibia  mass  to  be  more 
distributed  throughout  the  model.  It  can  be  noted  in  figure  3  that  using  the  midpoint  method  and 
an  increased  number  of  masses  leads  to  a  2.72%  decrease  in  the  maximum  tibia  force  as  well  as  a 
decrease  in  the  amplitude  of  oscillations. 


Figure  3.  Tibia  force  versus  time  for  the  27-spring,  midpoint  method  analytical  model. 

Using  an  implicit  numerical  method,  the  n-spring  model  allows  for  an  unlimited  distribution  of 
masses  and  number  of  springs  in  the  model.  The  use  of  the  two-step  Runge-Kutta  numerical 
integration  method  leads  to  a  further  decrease  of  0.34%  in  maximum  tibia  force  and  eliminates 
the  residual  oscillations  seen  in  both  the  forward  Euler  and  midpoint  method  solutions. 
Increasing  the  number  of  nodes  beyond  50  (figure  4)  proved  to  be  relatively  time-expensive  and 
unnecessary  since  the  solution  converged  with  this  many  nodes. 
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Figure  4.  Tibia  force  versus  time  for  the  50-spring,  two-step  Runge-Kutta  method  analytical  model. 

In  figure  5,  it  can  be  seen  that  as  the  number  of  springs  increases  in  the  midpoint  method  and  the 
two-step  Runge-Kutta  models,  the  calculated  maximum  tibia  force  converges.  However,  the 
calculated  force  is  approximately  twice  the  experimental  maximum  tibia  force  experienced 
without  boots. 
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Figure  5.  Force  versus  number  of  springs  for  different  numerical  methods. 

Note:  The  lines  show  minimum  and  maximum  no  boot  maximum  tibia  force  from  the  Bir  et  al.  (2006) 
measurements. 

In  figures  6  and  7,  the  calculated  force  maximums  from  DISP3  very  nearly  match  experimental 
data,  while  the  force  maximums  calculated  from  DISP2  demonstrate  that  the  floor  displacement 
data  resulted  from  a  much  more  severe  blast  than  experienced  in  the  TROSS  experiment. 
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Figure  6.  Scaled  maximum  tibia  force  versus  maximum  floor  deflection  . 


Note:  For  both  figures  6  and  7,  the  data  points  shown  use  different  numbers  of  springs  in  the  model;  the 
FE  model  uses  5  springs,  the  MM  model  uses  26  springs,  and  the  RK2  model  uses  50  springs. 
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Figure  7.  Scaled  maximum  tibia  force  versus  maximum  floor  velocity. 
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4.  GUI  Executable 


In  order  to  make  the  computation  more  user-friendly,  a  graphical  user  interface  (GUI)  executable 
was  developed  to  allow  the  user  to  reproduce  the  results  of  this  report  without  extensive 
knowledge  of  the  code.  A  screenshot  is  provided  in  figure  8,  showing  the  results  of  using  Bir  et 
al.  (2006)  Condition  1  TROSS  data. 


Figure  8.  GUI  executable  screenshot. 


To  improve  the  accuracy  of  the  predicted  tibia  force,  the  model  was  extended  further  to  include 
the  presence  of  the  ankle  and  foot,  and  the  boot.  The  parameters  listed  in  table  3,  derived  from 
the  Hybrid-Ill  LS-DYNA  model,  were  used  in  the  calculation  of  the  tibia  force. 
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Table  3.  Lengths  and  spring  constants  of  the  H-II1  tibia,  ankle/foot,  and  boot. 


Tibia 

Ankle/Foot 

Boot 

Length  (m) 

0.4523 

0.06119 

0.02967 

Spring  Constant  (N/m) 

26,954,124.92 

kfibia 

2,898,600 

The  values  for  the  spring  constants  of  the  ankle  and  foot  and  of  the  boot  are  approximated.  The 
assumption  that  the  spring  constant  of  the  ankle  was  equal  to  the  spring  constant  of  the  tibia  is 
nearly  arbitrary;  the  spring  constant  of  the  boot  assumes  a  hard  ethylene  propylene  diene 
monomer  (EPDM)  rubber  of  8.6  MPa  at  100%  elongation  and  a  cross-sectional  area  of  0.01  m  . 

The  mass  parameters — mfemur,  mtibia,  and  mfoot — are  shown  in  table  4. 

Table  4.  Masses  of  the  Hybrid-Ill  femur,  tibia,  and  foot. 


mfemur  (kg) 

4.16503 

m,ibia  (kg)  (includes  ankle) 

2.44792 

mf00t  (kg)  (includes  boot) 

1.38669 

With  this  further  modification,  the  model  places  a  mass  at  the  ankle  and  at  the  heel  and  provides 
a  different  spring  constant  for  both  the  ankle  and  the  sole  of  the  boot.  The  significantly  lower 
spring  constant  of  the  EPDM  rubber  in  the  sole  compared  to  the  bone  structures  greatly 
diminished  the  axial  tibia  force. 

For  force  verification  purposes,  an  effort  to  match  the  maximum  floor  displacements  and 
velocities  to  the  experimental  values,  the  TROSS  data  from  Bir  et  al.  (2006)  is  recreated 
manually  from  the  paper  through  measurements  and  extrapolations.  The  predicted  force  versus 
time  for  condition  1  is  shown  in  figure  9. 
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Tibia  Force  v.  Time 


Figure  9.  Calculated  tibia  force  versus  time  for  Bir  et  al.  (2006)  condition  1. 

Repeating  the  calculation  for  each  of  the  conditions  provided  in  Bir  et  al.  (2006)  resulted  in  a 
more  accurate  approximation  of  the  calculated  tibia  force  using  a  50-spring,  50-mass  model 
(table  5). 


Table  5.  Results  from  the  model  using  Bir  et  al.  (2006)  conditions  1  to  3  floor  displacement  inputs. 


Max  Floor 

Max  Floor 

Calculated  Max 

Experimental  Max 

Error  = 

Displacement 

Velocity 

Tibia  Force 

Tibia  Force 

(Calculated  -  Experimental)  / 

Condition 

(m) 

(m/s) 

(N) 

(N) 

Experimental 

1 

0.0125 

1.7 

4439 

3709 

0.20 

2 

0.0218 

3.3 

8148 

7000 

0.16 

3 

0.0269 

4.3 

10800 

9984 

0.08 

Since  the  model  includes  boots,  the  results  in  table  5  are  compared  to  the  experimental  results  for 
conditions  1  to  3.  From  the  results,  a  scaling  factor  (experimental/predicted)  is  calculated  for 
each  condition  and  a  quadratic  polynomial  is  used  to  relate  the  scaling  factor  to  its  corresponding 
maximum  floor  displacement  and  maximum  floor  velocity.  This  polynomial  should  be  further 
improved  through  the  verification  of  more  data  points. 

Since  the  scaling  factor  produced  using  the  maximum  floor  velocity  over-predicts  the  maximum 
tibia  force  as  seen  in  figure  10,  these  scaling  factors  are  used  in  the  model  to  provide  a 
conservative  estimate. 

It  can  be  seen  in  figures  10  and  1 1  that  the  model  with  modifications  to  account  for  the  properties 
of  the  ankle,  foot,  and  boot  results  in  a  more  accurate  prediction  of  the  axial  tibia  force.  Since  the 
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correlation  between  the  experimental  data  and  the  scaled  analytical  solution  is  stronger  using  the 
maximum  floor  velocity,  the  GUI  presents  the  scaled  analytical  solution  using  the  scaling  factor 
detennined  by  the  maximum  floor  velocity. 


Maximum  Tibia  Force  v.  Maximum  Floor  Deflection 


Figure  10.  Boot  maximum  tibia  force  versus  maximum  floor  displacement. 


Maximum  Tibia  Force  v.  Maximum  Floor  Velocity 


-B-  — Experimental 
—©■  —  Analytical 
— © — Scaled  Analytical 


Figure  1 1 .  Boot  maximum  tibia  force  versus  maximum  floor  velocity. 


5.  Summary  and  Conclusions 


The  one-dimensional  n-spring  analytical  model  given  simulated  floor  displacement  data 
produced  maximum  tibia  forces  within  the  same  magnitude  of  the  experimental  data.  Using  the 
forward  Euler,  midpoint  method,  and  two-step  Runge-Kutta  numerical  integration  methods,  the 
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developed  script  calculated  a  numerical  solution  for  the  tibia  force  given  the  floor  displacement 
over  time.  A  scaling  factor  dependent  on  the  maximum  floor  velocity  was  applied  to  obtain 
analytical  results  that  very  nearly  matched  the  experimental  results.  While  the  model  has  been 
validated  with  the  experimental  data  presented  in  Bir  et  al.  (2006),  further  verification  of  the 
model  with  a  larger  data  set  is  necessary. 
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Appendix  A.  Figures  Referenced  from  Bir  et  al.  (2006) 


Table  A-l.  Results  from  testing  conducted  using  TROSS  system  and  Hybrid-Ill  surrogate. 


Hybrid  III 

Plate  displacement 
(mm) 

Peak  plate  velocity 
(in/s) 

Tibia  Force  -  Z 
(N) 

C'ouditioii  1  no  boot 

12.3 

2.0 

5970 

C'ouditiou  1  boot 

12.6 

2.0 

3709 

C'ouditioii  2  uo  boot 

21.8 

3.4 

10740 

Condition  2  boot 

20.5 

3.4 

7000 

Condition  3  boot 

27.1 

4.5 

9984 

Figure  A-l.  Comparison  of  linear  impactor  data  with  TROSS  testing  using  Hybrid-Ill 
lower  limb  for  Condition  1. 
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Figure  A-2.  Comparison  of  linear  impactor  data  with  TROSS  testing  using  Hybrid-Ill 
lower  limb  for  Condition  2. 


Figure  A-3.  Comparison  of  linear  impactor  data  with  TROSS  testing  using  Hybrid-Ill 
lower  limb  for  Condition  3. 
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